* script to create distance bin plots
* October 2019


* Speed changes for "leavers" and "remainers" on SF--LALB
local ves "Container" // "Other Cargo" // 
local bins kmh_b
local eca_ys 2009 2011 // 2009
local yvar dist_b
local routes LALB_SFBay 
local dir_options Both // Both // NW SE 
local sbinds 0 1 
local preposts 1 // 0
local treat_opts 1 
local yscale "-50(50)150"


foreach prepost of local preposts {
	foreach treat_opt of local treat_opts {
		local treat_tag ""
		if `treat_opt'==0 {
			local treat_tag "_exempt"
		}
	
		foreach sb of local sbinds {  
			foreach eca_y of local eca_ys {
						
				capture drop sample
				if `eca_y' == 2009 {
					local cut eca1
					local run t_eca09 	
					local bw 150
					local donut 0
				}
				if `eca_y'==2011 {
					local cut eca2
					local run t_eca11
					local bw 150
					local donut 0
				}	
				gen sample = (`run'>=-`bw') & (`run'<=`bw') & treat==`treat_opt' 	& main_sample==1
					
				* in/out of sb channel post policy 
				capture drop case_samp
				gen case_samp = ( (`cut'==1) & SBind==`sb' ) | (`cut'==0)
				
				capture drop tmp_*
				* flagging and counting those in sample post policy
				* positive if in sample post cutoff
				gen case_samp_post = ( (`cut'==1) & SBind==`sb' ) 
				egen tmp_cnt_samp = total(case_samp_post) , by(route_id ves_id)
				drop case_samp_post
				
				* Limit Sample to vessels in sample that we observe pre
					* so keep only vessels we observe leaving post policy
					* or only vessels we observe staying post policy 

				gen tmp_pre = (`run'>=-`bw') & (`run'<-`donut') & tmp_cnt_samp>0
				egen tmp_cnt_pre = total(tmp_pre) , by(route_id ves_id)
				gen tmp_post = (`run'<=`bw') & (`run'>`donut') & tmp_cnt_samp>0 
				egen tmp_cnt_post = total(tmp_post) , by(route_id ves_id)

				gen tmp_pre_post = tmp_cnt_post>0 & tmp_cnt_pre>0 
				
				local prepost_opt "inlist(tmp_pre_post,1,0)" 
				if `prepost' == 1 {
					local prepost_opt "inlist(tmp_pre_post,1)"
				}
										
				* update sample
				gen sample2 = (sample==1) & `prepost_opt' & (case_samp==1)
				replace sample =  sample2
				drop sample2
						
				tab `cut' if sample
								
				foreach route of local routes {
					foreach dir_tag of local dir_options {
						
						if "`dir_tag'" == "Both" {
							local dir_cond "inlist(dirnw,0,1)"
						}
						if "`dir_tag'"=="NW" {
							local dir_cond "inlist(dirnw,1)"
						}
						if "`dir_tag'"=="SE" {
							local dir_cond "inlist(dirnw,0)"
						}
							
						reghdfe `yvar'  i.`bins' `cut'#i.`bins'  `cut'#i.`bins' ///  #c.`run'
						if route_name=="`route'" & `dir_cond' & vesseltype_regstr=="`ves'"  & sample  ///
							, absorb(ves_id#i.`bins') vce(cluster ves_id) coefleg
						est store speed_bin	
						parmest , saving(tmp, replace) label 
						
						*tabstat USflag dwt , by(`cut')

						reg `yvar' i.`bins' if e(sample) & `cut'==0 , noc
						est store base	
						parmest , saving(tmp_base, replace) label 
							
						* simple version
						/*
						coefplot (base , label("Baseline") keep( *.`bins' ) rename("^([0-9]+).`bins'$" = \1 , regex) recast(bar) noci )   ///
								 (speed_bin , label("Change") keep( 1.`cut'#*.`bins'  ) rename("^1.`cut'#([0-9]+).`bins'$" = \1 , regex)  )    ///
										, vert  nooff  ytitle("Distance Traveled (km)") xtitle("Speed Bin (km/h)") /// rename(^[0-9|0-9bn].eca_ind#([0-9]) = \1 , regex)
								name(diff_`dir_tag'_`route'_eca`eca_y',replace) ylab(,nogrid) 
						
						graph export "rd_plots\speedbins_`yvar'_`ves'_`eca_y'_`route'_SB`sb'_prepost`prepost'`treat_tag'.eps", replace
						*/			
								
						preserve 
						use tmp, replace		
						split parm , p("#")
						drop if parm2=="" // | parm3~="" gets rid of everything but interactions
						gen eca_ind = subinstr(subinstr(parm1,".`cut'","",.),"b","",.)
						gen kmh_b = subinstr(subinstr(subinstr(parm2,".kmh_b","",.),"b","",.),"o","",.)
						destring eca_ind kmh_b , replace
						save tmp_deltas , replace

						use tmp_base , replace
						gen kmh_b = subinstr(subinstr(subinstr(parm,".kmh_b","",.),"b","",.),"o","",.)
						destring kmh_b , replace
						rename estimate	pre_mean
						keep kmh_b pre_mean
						save tmp_base , replace

						use tmp_deltas , replace
						merge m:1 kmh_b using tmp_base
						drop _merge
						drop if eca_ind == 0
						drop parm*

						* generate counterfactual profile
						egen pre_sum = sum(pre_mean)
						egen d_sum = sum(estimate)

						gen pre_shr = pre_mean / pre_sum

						* using all distance after slowest 200km
						gen pre_cum = sum(pre_mean)
						gen	pre_cum_adj = max( 0 , pre_cum - 200 ) 
						gen pre_adj = pre_cum_adj - pre_cum_adj[_n-1]
						gen pre_adj_shr = pre_adj/(pre_sum - 200)

						gen d_prop = d_sum*pre_shr
						gen d_prop_adj = d_sum*pre_adj_shr

						gen kmh_b_mid = (kmh_b[_n+1] - kmh_b)/2  + kmh_b
						replace kmh_b_mid=kmh_b if kmh_b_mid==.

						* changes in fuel use
						gen f0 = kmh_b_mid^2 * pre_mean 
						gen f = kmh_b_mid^2 * estimate
						gen f_cf_prop = kmh_b_mid^2 * d_prop
						gen f_cf_prop_adj = kmh_b_mid^2 * d_prop_adj

						egen sum_f0 = sum(f0)
						egen sum_f = sum(f)
						egen sum_f_prop = sum(f_cf_prop)
						egen sum_f_prop_adj = sum(f_cf_prop_adj)

						gen ratio = sum_f / sum_f_prop_adj

						gen f_inc_per = sum_f / sum_f0
						gen f_inc_per_prop_adj = sum_f_prop_adj / sum_f0

						twoway  (bar pre_mean kmh_b , barw(2)) (bar d_prop_adj kmh_b , barw(2) bcolor(blue) ) /// 
							   || (rcap max95 min95 kmh_b , lcolor(dkorange)) || (scatter estimate kmh_b , mcolor(dkorange))  ///
							   , ytitle("Distance (km)") xtitle("Speed Bin (km/h)") ///
							   legend(cols(3) order(1 "Baseline" 4 "Change" 2 "CF Change" ) ) /// 
							   ylab(`yscale',nogrid) yline(0, lcolor(black))
						graph export "rd_plots\speedbins2_`yvar'_`ves'_`eca_y'_`route'_SB`sb'_prepost`prepost'.eps", replace
						
						restore 
						
					}
				}
			}
		}
	}
}

